.. _tutorial_1: TUTORIAL 1 : NVT Lennard-Jones fluid ==================================== :Authors: Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk), James Grant (r.j.grant{at}bath.ac.uk), and John Purton (john.purton{at}stfc.ac.uk) MC method in outline --------------------- In this tutorial we will perform the most straightforward Monte Carlo (MC) simulation - a simulation of a Lennard-Jones fluid in the canonical (NVT) ensemble. The triplet **NVT** signifies that number of particles, volume (hence, density) and temperature are kept constant during the entire length of this type of simulation. In this case only the particles are allowed to move, and a random walk is constructed following the approach of Metropolis *et al.* [#F1]_. The Metropolis algorithm iterates the following three steps that together constitute a single MC move (or attempt): 1. Select a particle at random and calculate its energy :math:U(\mathbf{r}_1). 2. Displace the particle randomly to a *trial* position and calculate its new energy :math:U(\mathbf{r}_2). 3. Accept or reject the particle displacement according to the Metropolis acceptance rule: .. math:: P_{\mathrm{acc}}(\mathbf{r}_1 \rightarrow \mathbf{r}_2) = \min(1, \exp \{- \beta [U(\mathbf{r}_2) - U(\mathbf{r}_1)] \} ) How do we achieve this? ----------------------- **Ensure detailed balance**: it is crucial to ensure a truly random (uniform) selection of trial particle positions! Otherwise an unwanted, hidden bias can be introduced into the MC sampling procedure, which would invalidate the obtained statistics. The common way to move (translate) a particle is to generate a random displacement vector :math:\Delta\mathbf{r} bound within a cube with a given side legnth :math:\Delta_{max} (the so-called *displacement parameter*) and add it to the current position of the particle under trial. .. sphere of a given radius :math:\Delta_{max} Then, pick up a random number between 0 and 1 and accept or reject the move as illustrated below: .. figure:: ./images/Metropolis-acceptance.jpg :width: 500px DL_MONTE Input Files -------------------- The purpose of this tutorial is to introduce the basic file formats that are used by DL_MONTE and then run an example MC simulation. DL_MONTE requires at least 3 standard input files (similar to DL_POLY): FIELD -- units, atoms, molecules, topology and force-field specs (fixed format) CONFIG -- initial configuration for the system, including cell vectors (fixed format) CONTROL -- simulation parameters: external conditions, options and flags (free format) NOTE: The file names are *fixed*, i.e. only these file names are recognised by DL_MONTE. In some more intricate cases other input files might be needed (e.g. for restarting a simulation for continuation). We will introduce these files and discuss each of them in turn. Full details of their structure, formats and default values can be found in the DL_MONTE manual. To check the input files for the current tutorial, navigate to the corresponding directory:: [DL_MONTE-2_tutorial]$cd exercises/tutorial_1 [DL_MONTE-2_tutorial]$ ls -l NOTE: The first line in all the three files is an arbitrary *title* which can contain up to 80 characters (any symbol over that limit will be ignored by DL_MONTE). The titles can be different in each file. FIELD ----- In DL_MONTE the system is abstracted into atoms and molecules: atoms are treated as point-like particles, and molecules are collections of atoms which can be moved collectively. This, along with a versatile selection of potential forms which can be combined into different force fields, allows for simulation of a wide range of systems comprised of possible combinations of free unconnected atoms (so-called atomic/ionic fields), and molecules possessing structure: both rigid and flexible. The FIELD file defines all the atomic and molecular types present in the system, but also possibly the molecular structure (so-called topology: bonds, angles, etc). Finally, it also specifies the interatomic and (optional) external potentials. A detailed description of all the available force-field parameters and their interplay can be found in the DL_MONTE manual. The file structure is similar to that used by DL_POLY, but there are some important differences owing to the requirements for Grand-Canonical simulations (considered later). The FIELD file for our first tutorial is shown below .. code-block:: html :linenos: Lennard-Jones, 2.5*sigma cut-off, sigma = 1 Angstrom, epsilon = 1 internal unit = 10 J/N_A CUTOFF 2.5 # cut-off radius (must be the 2-nd line - after the title) UNITS internal # energy unit (must be the 3-rd line - after CUTOFF) NCONFIGS 1 # number of configurations in CONFIG file (must be the 4-th line) ATOM TYPES 1 # number of atom types or elements (must be the 5-th line) LJ core 1.0 0.0 # mass and charge of the atom(s) MOLECULE TYPES 1 # number of molecular types lj # molecule name (case sensitive) MAXATOM 512 # max number of atoms within current molecular type FINISH # closing the molecular type specs VDW 1 # number of short-ranged Van der Waals potentials LJ core LJ core lj 1.0 1.0 # atom names and types for a VDW pair, LJ epsilon and sigma CLOSE # closing FIELD The important points to note are: the *cutoff = 2.5* for the short-ranged VDW interactions is in Angstroms the energy *unit* used in this case is *internal = 10 J/mol* [or 0.01 kJ/mol], the CONFIG file contains *NCONFIGS = 1* configuration. the only *ATOM* type present has the name *LJ* and the type *core*, having *mass = 1.0* and *charge = 0.0*. All the atoms in the system are contained within a single molecule called *lj*. The number of atoms in this molecule is restricted to 512 at most. Note that the CONFIG file cannot contain more than that number of atoms within that type of a molecule. As a matter of fact, for an NVT simulation where the number of particles cannot change, the CONFIG file *must* contain the total number of atoms that is specified in the FIELD file! The atoms interact through a Lennard-Jones interaction potential, which is specified by the keyword *lj*, and the LJ parameters are set to: *\epsilon = 1.0* and *\sigma = 1.0*. .. Unlike DL_POLY there is no one-to-one correspondence between the FIELD and CONFIG file thus we only need to specify the particle type under the keyword ATOMS. CONFIG ------ The CONFIG file contains the starting configuration specification. It defines the cell type and geometry, sets the cell vectors (in the matrix form), the actual and maximum numbers of molecules, atoms within each molecule, and also the coordinates of all the atoms. Below we illustrate the structure of the CONFIG file with a sample of the configuration provided for this tutorial. .. code-block:: html :linenos: Exercise 1: LJ fluid in NVT ensemble 0 0 11.7452 0.00000 0.00000 0.00000 11.7452 0.00000 0.00000 0.00000 11.7452 NUMMOL 1 1 MOLECULE lj 512 512 LJ core 0 0 0 LJ core 0 0 0.125 LJ core 0 0 0.25 LJ core 0 0 0.375 ...... In line 2 the first number is an integer switch operating in the same way as *levcfg* in DL_POLY. It should be set to zero if the CONFIG file only contains the particle coordinates (which is sufficient for starting an MC simulation). However, if non-zero, this flag allows DL_MONTE to read in CONFIG files originating from MD (DL_POLY) simulations and containing (x,y,z) components of the particle velocities (*levcfg = 1*) or velocities and forces (*levcfg = 2*) which would appear on the lines directly following the particle coordinates. Note, however, that DL_POLY CONFIG files cannot be used directly, without certain modifications, as input for DL_MONTE simulations. The second integer flag determines the convention used for the particle coordinates: *0 for fractional* representaion, and *non-zero for Cartesian*: 1 == cubic, 2 == orthorhombic, 3 == palleliped. These are followed by the cell matrix. In our case the CONFIG file contains fractional coordinates and, since the cell matrix is diagonal and all its diagonal elements are equal, the simulation cell is cubic. The line starting with *NUMMOL* keyword specifies the *total number of molecules present*, followed by the *maximum number of molecules allowed for each molecular type* defined in the FIELD file. In this case we have only one molecular type *lj*, and there is only one molecule of this type in the CONFIG file. The molecule comprises 512 atoms (LJ particles) which equals the maximum number. Finally, the atoms' specs are read in. The atom data follows and for each particle, its name (usually chemical symbol) and type (core, semi and metal that can be abbreviate to c, s, and m) are specified. On the next line these are followed by the coordinates in a format consistent with that given on line 2. .. figure:: ./images/tutorial1-CONFIG1.png :width: 450px CONTROL ------- The CONTROL provides directives to DL_MONTE how to undertake the calculations and switches on or off functionality. The CONTROL file in this example is: .. code-block:: html :linenos: NVT simulation of Lennard-Jones fluid #use ortho finish seeds 12 34 56 78 # Seed RNG seeds explicitly to the default nbrlist auto # Use a neighbour list to speed up energy calculations maxnonbondnbrs 512 # Maximum number of neighbours in neighbour list temperature 1.4283461511745 # T* = 1.1876 = T(K) * BOLTZMAN (see constants_module.f90) steps 10000 # Number of moves to perform in simulation equilibration 0 # Equilibration period: statistics are gathered after this period print 1000 # Print statistics every 'print' moves stack 1000 # Size of blocks for block averaging to obtain statistics sample coord 10000 # How often to print configurations to ARCHIVE.000 revconformat dlmonte # REVCON file is in DL_MONTE CONFIG format archiveformat dlpoly4 # ARCHIVE.000/HISTORY.000/TRAJECTORY.000 format # In this case: HISTORY.000 in DLPOLY4 style move atom 1 100 # Move atoms with a weight of 100 LJ core start The first line is the title and the second contains the keyword *finish*. We will see later that there a number of directives in DL_MONTE that can only occur before this directive which *must* be present in the CONTROL file. *Seeds* followed by a series of 4 integers provides a reproducible seed, otheriwse *ranseed* generates a random seed from the system clock at initialisation. The diretives *nbrlist* and *maxnonbondnbrs* control the size and administration of the neighbourlist used by DL_MONTE to optimise performance, and are explained further in one of the exercises. *Temperature* (K) is self explanatory while *steps* is the length of the simulation in attempted moves. NOTE: Providing that LJ epsilon = 1.0 in FIELD, the reduced LJ temperature T* = T(K)*R, where R = BOLTZMAN (see constants_module.f90 or beta = 1/RT(K) in OUTPUT.000). In general, T* = RT(K)/epsilon = 1/(epsilon*beta), where epsilon is in *internal energy units*. *Equilibration* etc control the detail of how data is output from DL_MONTE. The *revconformat dlmonte* instruction describes the format of the output file REVCON.000 which contains the final configuration of the simulation and can be used for continuing a simulation, *dlpoly2*, *dlpoly4* are other options. *archiveformat dlpoly4* describes the format of the trajectory file here it will be HISTORY.000, equivalent of the HISTORY in DLPOLY4. The *dlpoly{2,4}* formats are readable by common visualisation packages such as vmd. Full options are detailed in the manual. The directive *move atoms* is where we begin to touch on the fundamental control of the simulation. The key feature here is that DL_MONTE will not do anything unless told to do so (N.B. While this gives DL_MONTE great flexibility it means also means that it may be possible to ask DL_MONTE to perform ill-defined calculations). In this simple *NVT* ensemble only the particles move, thus *move atoms* tells DL_MONTE to move *1* atom type, with a weight of *100*. The line(s) following this detail the atom and type. Finally the *start* directive ends the *CONTROL* file and instructs DL_MONTE to start the simulation. We are now ready to run DL_MONTE:: [tutorial_1]$DL_MONTE.SRL.X & When running DL_MONTE locally (i.e. in the "interactive" mode), to check if the simulation is still running, we can use the **top** command:: [tutorial_1]$ top which produces a list of the active processes. If DL_MONTE.SRL.X is still running its name will be found in the list. Press **q** key to quit the interactive top output. Examining the output files -------------------------- A successful DL_MONTE calculation will produce a number of output files:: * OUTPUT.000 - details of the simulation, statistics, running time, or errors if the calculation failed. * REVCON.000 - the final configuration in the format specified * PTFILE.000 - statistics though will eventually be deprecated in favour of * YAMLDAT.000 - statistics in the yaml format * ARCHIVE.000/HISTORY.000/TRAJECTORY.000 - the trajectory in the specified format For analysis we will typically process the YAMLDAT.000 and visualise the trajectory files. However for understanding how the simulation proceeds it is useful to have some familiarity with the OUTPUT file. The file begins with a header detailing the version, authors and suggested citations, followed by the brief summary of details of the simulation as specified in the input files. A section headed *simulation parameters* then specifies all parameter values that will be used within the simulation. The final step before starting the calculation is to determine the initial energy of the system and the details of this are printed in a block:: -------------------------------------------------- initial energies -------------------------------------------------- break down of energies for box: 1 total energy -0.7037212307E+03 reciprocal space coulomb 0.0000000000E+00 real space coulomb 0.0000000000E+00 external mfa coulomb 0.0000000000E+00 nonbonded two body (vdw) -0.7037212307E+03 bonded two body (pair) 0.0000000000E+00 nonbonded three body 0.0000000000E+00 bonded three body (angle) 0.0000000000E+00 bonded four body (angle) 0.0000000000E+00 many body energy 0.0000000000E+00 external potential energy 0.0000000000E+00 total virial 0.0000000000E+00 volume 0.1620247087E+04 This is followed by a partial breakdown per molecule type in the system and the time taken to initialise the calculation. There after every *print* steps as specified in the CONTROL file a block is printed to the OUTPUT.000 file:: Iteration 10000 - elapsed time: 0 h : 0 m : 0.301 s ==================================================================================================== step en-total h-total coul-rcp coul-real step en-vdw en-three en-pair en-angle step en-four en-many en-external en-extMFA step volume cell-a cell-b cell-c step alpha beta gamma r-av en-total h-total coul-rcp coul-real r-av en-vdw en-three en-pair en-angle r-av en-four en-many en-external en-extMFA r-av volume cell-a cell-b cell-c r-av alpha beta gamma ---------------------------------------------------------------------------------------------------- 10000 -0.1006733261E+04 -0.3030120301E+03 0.0000000000E+00 0.0000000000E+00 -0.1006733261E+04 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.1620247087E+04 0.1174520000E+02 0.1174520000E+02 0.1174520000E+02 0.9000000000E+02 0.9000000000E+02 0.9000000000E+02 -0.1014992775E+04 -0.3112715438E+03 0.0000000000E+00 0.0000000000E+00 -0.1014992775E+04 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.1620247087E+04 0.1174520000E+02 0.1174520000E+02 0.1174520000E+02 0.9000000000E+02 0.9000000000E+02 0.9000000000E+02 LJ c 512.0000 512.0000 lj 1.0000 1.0000 ---------------------------------------------------------------------------------------------------- equilibration period ended at step 10000 ---------------------------------------------------------------------------------------------------- This specifies the same breakdown in a tabulated form, including the iteration number and time taken. By using the command:: [tutorial_1]$grep time OUTPUT.000 you can quickly gauge the progress of your calculation. DL_MONTE runs a *check* of the system at regular intervals to verify that the system is behaving correctly. This typically involves calculating the energy from scratch and comaring it with the running total, the result is printed to the OUTPUT.000 with a line looking like:: [tutorial_1]$ grep U_recalc OUTPUT.000 Workgroup 0, box 1 check: U_recalc - U_accum = 0.35129E-10 0.29580E-10 0.57773E-13 -0.34894E-13 (internal, kT, kT/atom, dU/U) The final value is the relative energy difference which should be of the order of the working precision, typically ending E-13 or E-14. Finally at the end of the simulation a summary block detailing the average values during the simulation (excluding the first *equilbration* steps) and their fluctations, the final energies, and 'processing data' detailing time, move data and final parameters:: ---------------------------------------------------------------------------- processing data ---------------------------------------------------------------------------- total no of attempted & empty atom moves : 110000 0 successful no of atom moves & acceptance : 42987 0.39079091 displacement (Angstroms) for LJ c : 0.4160 pure simulation time (excl. init.): 0 h : 0 m : 3.492 s total elapsed time (incl. init.): 0 h : 0 m : 3.498 s normal exit Upon a successful completion of a simulation run, the OUTPUT.000 file will end with the string **normal exit**. To quickly check for this try:: [tutorial_1]$grep "normal exit" OUTPUT.000 Basic analysis of the output ---------------------------- .. To see if the simulation converged to equilibrium, add the following line:: .. yamldata 1000 .. to the CONTROL file and **rerun the simulation**. To examine the energy evolution along the MC time, first run the script:: [tutorial_1]$ strip_yaml.sh energy and then plot the energy vs MC time:: [tutorial_1]$gnuplot gnuplot> plot './energy.dat' u 1:2 w l t "Energy(MC step)" .. figure:: ./images/tutorial1-energy.jpg :width: 640px To visualise the trajectory we can use vmd:: [tutorial_1]$ vmd To visualise the trajectory:: File New molecule Determine file type -> Browse and Select HISTORY.000 .. figure:: ./images/tutorial1-VMD.png :width: 900px Select the correct format:: -> DLPOLY V3 History Note that in a newer version of VMD it is called DLPOLY-4 .. figure:: ./images/tutorial1-VMD2.png :width: 400px Questions to ask yourself: ^^^^^^^^^^^^^^^^^^^^^^^^^^ * What was the initial configuration, and what did happen during simulation? * Did the system reach equilibrium? * Can we improve on the stats? How? Extra exercises to try and do ----------------------------- Now try the extension exercises to learn more about funcitonality within DL_MONTE and to optimise your calculation: :ref:tut1_ex1 :ref:tut1_ex2 :ref:tut1_ex3 With each exercise we recommend you create a copy of the inputs in a sub-directory:: [tutorial_1]$mkdir ex1 [tutorial_1]$ cp CONFIG CONTROL FIELD ex1 [tutorial_1]\$ cd ex1 .. Link to next tutorial Or move on to :ref:tutorial_2 and NPT simulations. .. rubric:: Footnotes .. [#f1] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.N. Teller, and E. Teller. Equation of state calculations by fast computing machines. J. Chem. Phys. , 21:1087 1092, 1953.